rEDM should be the version 1.2.3
https://cran.r-project.org/src/contrib/Archive/rEDM/
library(mathjaxr)
library(rEDM); packageVersion("rEDM")
We are currently building a new version of rEDM, and a version for python pyEDM. Be on the look out for version 1.0 for rEDM. Updates can be found at https://github.com/SugiharaLab
[1] ‘1.2.3’
library(parallel)
library(doParallel)
Loading required package: foreach
Loading required package: iterators
library(foreach)
library(Kendall)
library(MASS)
library(dplyr)
Attaching package: ‘dplyr’
The following object is masked from ‘package:MASS’:
select
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(glmnet)
Loading required package: Matrix
Loaded glmnet 4.1-8
library(ggplot2)
#source("../Demo_MDR_function.R")
source("./Demo_MDR_function2.R") #from 2024/05/07
4th-order explicit Runge-Kutta method with a fixed interval function
rk4 <- function(in_vec, ref_data = NULL, time, h_interval, dim, diff_vec)
{
#in_vec: the vector with the value at "time"
#out_vec: the vector with which the updated vector value after one time step is saved
#time: time value
#h_interval: time step for discritizing ODE
#dim: dimension of the ODE
#diff_vec: function that determines the r.h.s. of the ODE
k1 <- numeric(dim)
k2 <- numeric(dim)
k3 <- numeric(dim)
k4 <- numeric(dim)
temp_vec <- numeric(dim)
h_half <- h_interval / 2.0
t_h <- time + h_half
k1 <- diff_vec(in_vec, ref_data, time, h_interval, dim) # calculate k1
temp_vec <- in_vec + h_half * k1
k2 <- diff_vec(temp_vec, ref_data, t_h, h_interval, dim) # calculate k2
temp_vec <- in_vec + h_half * k2
k3 <- diff_vec(temp_vec, ref_data, t_h, h_interval, dim) # calculate k3
temp_vec <- in_vec + h_interval * k3
k4 <- diff_vec(temp_vec, ref_data, time + h_interval, h_interval, dim) # calculate k4
temp_vec <- in_vec + (h_interval / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
temp_vec
}
Interaction strength estimation for time-discrete model (DE) framework
2-Host 2-Parasitoid System with natural mortality (i.e., with overlapping generations)
Note: When the natural mortality goes to infinity, the model converges to 2-Host 2-Parasitoid Nicolson-Baily Model
\[ H_{k,t+1} = H_{k,t} + r_k exp[-a_{k}H_{k,t} - c_{1k} P_{1,t} - c_{2k} P_{2, t}] H_{k,t} - (1 - exp[-a_{k}H_{k,t} - c_{1k} P_{1,t} - c_{2k} P_{2,t} - m_{H,k}]) H_{k,t} \quad (k=1,2) \tag{D1} \] It is set as \[ H_{k,t+1} = H_{k,t} + F_{Hk}(H_{1,t}, H_{2,t}, P_{1,t}, P_{2,t}) \]
\[ P_{j,t+1} = P_{j,t} + \sum_{k=1,2} exp(-a_k H_{k,t}) (1 - exp[-c_{1k} P_{1,t} - c_{2k} P_{2,t}]) \frac{c_{jk} P_{j,t}}{c_{1k} P_{1,t} + c_{2k} P_{2,t}} H_{k,t} -(1 - exp(-m_{P,j})) P_{j,t} \quad (j=1,2) \tag{D2} \] It is set as \[ P_{j,t+1} = P_{j,t} + F_{Pj}(H_{1,t}, H_{2,t}, P_{1,t}, P_{2,t}) \]
r1 <- 2.0
r2 <- 2.0
a1 <- 0.1
a2 <- 0.1
c11 <- 0.3 #P1 -> H1
c21 <- 0.1 #P2 -> H1
c12 <- 0.1 #P1 -> H2
c22 <- 0.3 #P2 -> H2
mH1 <- 0.1
mH2 <- 0.1
mP1 <- 0.1
mP2 <- 0.1
j_Host1 <- 1
j_Host2 <- 2
j_Paras1 <- 3
j_Paras2 <- 4
Note that FH1, FH2, FP1, FP2, are defined as the changes of the abundance (X_t+1 - X_t)
FH1 <- function(in_vec, t) {
term1 <- r1*exp(-a1 * in_vec[j_Host1] - c11 * in_vec[j_Paras1] - c21 * in_vec[j_Paras2]) * in_vec[j_Host1]
term2 <- (1.0 - exp(-a1 * in_vec[j_Host1] - c11 * in_vec[j_Paras1] - c21 * in_vec[j_Paras2] - mH1)) * in_vec[j_Host1]
return(term1 - term2)
}
FH2 <- function(in_vec, t) {
term1 <- r2*exp(-a2 * in_vec[j_Host2] - c12 * in_vec[j_Paras1] - c22 * in_vec[j_Paras2]) * in_vec[j_Host2]
term2 <- (1.0 - exp(-a2 * in_vec[j_Host2] - c12 * in_vec[j_Paras1] - c22 * in_vec[j_Paras2] - mH2)) * in_vec[j_Host2]
return(term1 - term2)
}
FP1 <- function(in_vec, t) {
term1 <- exp(-a1 * in_vec[j_Host1]) * (1 - exp(-c11 * in_vec[j_Paras1] - c21 * in_vec[j_Paras2])) * c11 * in_vec[j_Paras1] / (c11 * in_vec[j_Paras1] + c21 * in_vec[j_Paras2]) * in_vec[j_Host1]
term2 <- exp(-a2 * in_vec[j_Host2]) * (1 - exp(-c12 * in_vec[j_Paras1] - c22 * in_vec[j_Paras2])) * c12 * in_vec[j_Paras1] / (c12 * in_vec[j_Paras1] + c22 * in_vec[j_Paras2]) * in_vec[j_Host2]
term3 <- (1.0 - exp(- mP1)) * in_vec[j_Paras1]
return(term1 + term2 - term3)
}
FP2 <- function(in_vec, t) {
term1 <- exp(-a1 * in_vec[j_Host1]) * (1 - exp(-c11 * in_vec[j_Paras1] - c21 * in_vec[j_Paras2])) * c21 * in_vec[j_Paras2] / (c11 * in_vec[j_Paras1] + c21 * in_vec[j_Paras2]) * in_vec[j_Host1]
term2 <- exp(-a2 * in_vec[j_Host2]) * (1 - exp(-c12 * in_vec[j_Paras1] - c22 * in_vec[j_Paras2])) * c22 * in_vec[j_Paras2] / (c12 * in_vec[j_Paras1] + c22 * in_vec[j_Paras2]) * in_vec[j_Host2]
term3 <- (1.0 - exp(- mP2)) * in_vec[j_Paras2]
return(term1 + term2 - term3)
}
#function to calculate all changes four dimensional
diff_4HP <- function(in_vec, t, dim){
temp_vec <- numeric(dim)
temp_vec[j_Host1] <- FH1(in_vec, t)
temp_vec[j_Host2] <- FH2(in_vec, t)
temp_vec[j_Paras1] <- FP1(in_vec, t)
temp_vec[j_Paras2] <- FP2(in_vec, t)
temp_vec
}
dim_model2 <- 4
nv_4HP0 <- numeric(dim_model2)
nv_4HP <- numeric(dim_model2)
nv_4HP0[j_Host1] <- 1.0
nv_4HP0[j_Host2] <- 1.0
nv_4HP0[j_Paras1] <- 0.1
nv_4HP0[j_Paras2] <- 0.1
t <- 0 # initial condition (initial time, 0)
end_time <- 1000
#Setting initial condition
nv_4HP[j_Host1] <- 1.0
nv_4HP[j_Host2] <- 0.0
nv_4HP[j_Paras1] <- 0.1
nv_4HP[j_Paras2] <- 0.0
# write initial condition
cat(t, nv_4HP, "\n")
0 1 0 0.1 0
# record the initial condition
HP_result <- data.frame(t = t, H1 = nv_4HP[j_Host1], H2 = nv_4HP[j_Host2], P1 = nv_4HP[j_Paras1], P2 = nv_4HP[j_Paras2])
#For the transient dynamics (RUN BURN)
for(i in 1:end_time){
nv_4HP <- nv_4HP + diff_4HP(nv_4HP, t, dim_model2)
t <- t + 1 #update time
HP_result <- rbind.data.frame(HP_result, c(t, nv_4HP[j_Host1], nv_4HP[j_Host2], nv_4HP[j_Paras1], nv_4HP[j_Paras2]))
}
head(HP_result)
#Trial of the plot
plot(HP_result$t, HP_result$H1, type = "l", ylim = c(0,20))
par(new = T)
plot(HP_result$t, HP_result$P1, type = "l", ylim = c(0,20), col = "red")
t <- 0 # initial condition (initial time, 0)
end_time <- 2000
#Setting initial condition; it is important to start with asymmetric densities to reach quickly the chaotic attractor
nv_4HP[j_Host1] <- 1.0
nv_4HP[j_Host2] <- 0.5
nv_4HP[j_Paras1] <- 0.1
nv_4HP[j_Paras2] <- 0.2
# write initial condition
cat(t, nv_4HP, "\n")
0 1 0.5 0.1 0.2
# record the initial condition
HP_result <- data.frame(t = t, H1 = nv_4HP[j_Host1], H2 = nv_4HP[j_Host2], P1 = nv_4HP[j_Paras1], P2 = nv_4HP[j_Paras2])
#For the transient dynamics (RUN BURN)
for(i in 1:end_time){
nv_4HP <- nv_4HP + diff_4HP(nv_4HP, t, dim_model2)
t <- t + 1 #update time
HP_result <- rbind.data.frame(HP_result, c(t, nv_4HP[j_Host1], nv_4HP[j_Host2], nv_4HP[j_Paras1], nv_4HP[j_Paras2]))
}
head(HP_result)
#Trial of the plot
plot(HP_result$t, HP_result$H1, type = "l", ylim = c(0,5))
par(new = T)
plot(HP_result$t, HP_result$P1, type = "l", ylim = c(0,5), col = "red")
par(new = T)
plot(HP_result$t, HP_result$H2, type = "l", ylim = c(0,5), lty = "dashed")
par(new = T)
plot(HP_result$t, HP_result$P2, type = "l", ylim = c(0,5), col = "red", lty = "dashed")
DAY = Sys.Date()
saveRDS(HP_result, paste("HP4p_result_", DAY, ".obj", sep = ""))
Based on https://ushio-ecology-blog.blogspot.com/2019/12/20191225blogger0007.html
Loading the default data at 2024-07-25 and taking the final 200 data points in the time series
HP_result <- readRDS("HP4p_result_2024-07-25.obj")
HP_result <- HP_result[1801:2000, ]
HP_result.mean <- apply(HP_result[, -1], 2, mean, na.rm = T) # mean abundance
HP_result.sd <- apply(HP_result[, -1], 2, sd, na.rm = T) # SD of abundance
mean_mat <- matrix(HP_result.mean, nrow = nrow(HP_result) , ncol = length(HP_result.mean), byrow = TRUE)
sd_mat <- matrix(HP_result.sd, nrow = nrow(HP_result) , ncol = length(HP_result.mean), byrow = TRUE)
HP_result_s <- (HP_result[, -1] - mean_mat) /sd_mat # Standardized data set
head(HP_result_s)
# data block for smap
HP_smap_block_s <- data.frame(H1 = HP_result_s$H1, H2 = HP_result_s$H2, P1 = HP_result_s$P1, P2 = HP_result_s$P2)
# delete a few column
# conducting univariate smap for obtaining the optimal theta
HP_m_smap_res_s <- block_lnlp(HP_smap_block_s, method = "s-map",
target_column = j_Host1,
theta = c(0, 1e-04, 3e-04, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8),
silent = TRUE)
Plotting the effect of theta
plot(HP_m_smap_res_s$theta, HP_m_smap_res_s$rmse, type = "b")
(HP_opt_theta_s <- HP_m_smap_res_s$theta[which.min(HP_m_smap_res_s$rmse)])
[1] 8
Conducting multivariate smap with the optimal theta
HP_m_smap_res2_s <- block_lnlp(HP_smap_block_s, method = "s-map",
target_column = j_Host1,
theta = HP_opt_theta_s,
silent = TRUE,
save_smap_coefficients = TRUE)
head(as.data.frame(HP_m_smap_res2_s$smap_coefficients[[1]]))
final_ID <- 200
delta_11 <- 1
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11, type = "l", col = 1, ylim = c(-0.9, 1.0), xlab = "time", ylab = "Smap coefficient: X -> Host1", main = "S-map for H1")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_2[-final_ID], type = "l", col = "red", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID], type = "l", col = "blue", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_4[-final_ID], type = "l", col = "green", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
legend("topright",
legend = c("S-map from H1 to H1 - delta11", "S-map from H2 to H1", "S-map from P1 to H1", "S-map from P2 to H1"),
col = c("black", "red", "blue", "green"),
lty = c(1, 1, 1, 1),
cex = 0.8)
gg_HP_result00 <- data.frame(time = HP_result$t[-final_ID], smap = HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID], type = rep("H1toH1", 199))
gg_HP_result01 <- data.frame(time = HP_result$t[-final_ID], smap = HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11, type = rep("H1toH1_adjusted", 199))
gg_HP_result02 <- data.frame(time = HP_result$t[-final_ID], smap = HP_m_smap_res2_s$smap_coefficients[[1]]$c_2[-final_ID], type = rep("H2toH1", 199))
gg_HP_result03 <- data.frame(time = HP_result$t[-final_ID], smap = HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID], type = rep("P1toH1", 199))
gg_HP_result04 <- data.frame(time = HP_result$t[-final_ID], smap = HP_m_smap_res2_s$smap_coefficients[[1]]$c_4[-final_ID], type = rep("P2toH1", 199))
gg_HP_result <- rbind.data.frame(gg_HP_result00, gg_HP_result01)
gg_HP_result <- rbind.data.frame(gg_HP_result, gg_HP_result02)
gg_HP_result <- rbind.data.frame(gg_HP_result, gg_HP_result03)
gg_HP_result <- rbind.data.frame(gg_HP_result, gg_HP_result04)
fig1_ggplot <- ggplot(data = gg_HP_result) +
scale_colour_manual(values = c("black", "#ff4b00","#4dc4ff", "#f6aa00", "#804000")) +
scale_shape_manual(values = c(1, 15, 16, 17, 5)) +
 geom_line(aes(x = time, y = smap, group = type, color = type)) +
geom_point(aes(x = time, y = smap, color = type, shape = type), size = 2) +
xlim(1800, 1999) + ylim(-0.5, 1.5) + theme_bw() + theme(legend.position = "bottom") +
labs(title = "Figure 1: Time evolution of S-map coefficients", x = "time", y = "Unadjusted and Adjusted S-map coefficients")
plot(fig1_ggplot)
Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).
ggsave("fig1_ggplot.pdf", width = 10, height = 5)
Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).
Applying MDR S-map to the 4sp coupled host-parasitoid model
Read dataset
da.range <- 1801:2000 # Subsample for data analysis
out.sample <- T # T/F for out-of-sample forecast
nout <- 2 # number of out-of-sample
da.name <- '4p_model'
do <- readRDS("HP4p_result_2024-07-25.obj")
dot <- do[da.range, 1] # data time
do <- do[da.range, -1] # subset of data without t
ndo <- nrow(do)
nin <- ndo - nout # library sample size
Standardization
# In-sample
do.mean <- apply(do[1:nin, ], 2, mean, na.rm = T) # mean abundance in in-sample
do.sd <- apply(do[1:nin, ], 2, sd, na.rm = T) # SD of abundance in in-sample
d <- do[1:(nin-1), ] # In-sample dataset at time t
d_tp1 <- do[2:(nin), ] # In-sample dataset at time t+1
ds <- (d - repmat(do.mean, nrow(d), 1))*repmat(do.sd, nrow(d), 1) ^ -1 # Normalized in-sample dataset at time t
ds_tp1 <- (d_tp1 - repmat(do.mean, nrow(d_tp1), 1))*repmat(do.sd, nrow(d_tp1), 1) ^ -1 # Normalized in-sample dataset at time t+1
# Out-sample
if(out.sample & nout != 0){
d.test <- do[nin:(ndo - 1), ] # Out-of-sample dataset at time t
dt_tp1 <- do[(nin + 1): ndo, ] # Out-of-sample dataset at time t+1
ds.test <- (d.test - repmat(do.mean, nrow(d.test), 1))*repmat(do.sd, nrow(d.test), 1) ^ -1 # Normalized out-of-sample dataset at time t
dst_tp1 <- (dt_tp1 - repmat(do.mean, nrow(dt_tp1), 1))*repmat(do.sd, nrow(dt_tp1), 1) ^ -1 # Normalized out-of-sample dataset at time t+1
}else{d.test <- dt_tp1 <- ds.test <- NULL}
# Compiled data at time t
ds.all <- rbind(ds, ds.test)
head(ds.all)
Initialization of pseudo-random number generator
seed <- 49563
set.seed(seed)
Find the optimal embedding dimension, based on univariate simplex projection
Emax <- 4 # equal to or smaller than the number of nodes (dimension of the DE system)
cri <- 'rmse' # model selection
Ed <- NULL
forecast_skill_simplex <- NULL
for(i in 1:ncol(ds)){
spx.i <- simplex(ds[, i], E = 2:Emax)
Ed <- c(Ed, spx.i[which.min(spx.i[, cri])[1], 'E'])
forecast_skill_simplex <- c(forecast_skill_simplex, spx.i[which.min(spx.i[, cri])[1], 'rho'])
}
Ed # The optimal embedding dimension for each variable
[1] 2 3 2 2
Optimal embedding dimension and forecasting skill
Ed
[1] 2 3 2 2
forecast_skill_simplex
[1] 0.9882117 0.9882610 0.9917444 0.9875901
Find causal variables by CCM analysis that will be used for multiview embedding.
Warning: It is time consuming for calculating the causation for each node
CCM causality test for all node pairs (the function was updated at 2024/05/07)
ccm.out <- ccm.fast.demo(ds, Epair = T, cri = cri, Emax = Emax)
variable 1 ccm completed: 4.148 sec
variable 2 ccm completed: 4.128 sec
variable 3 ccm completed: 4.369 sec
variable 4 ccm completed: 4.172 sec
ccm.sig <- ccm.out[['ccm.sig']]
ccm.rho <- ccm.out[['ccm.rho']]
#To avoid overwrite the original files, we save them with different names,identified by date.
DAY = Sys.Date()
saveRDS(ccm.sig, paste('ccm_sig', da.name,'nin', nin, DAY, '.obj', sep = '_'))
saveRDS(ccm.rho, paste('ccm_rho', da.name,'nin', nin, DAY, 'demo.obj', sep = '_'))
Just in case not conducting CCM but just loading the result of CCM at 2024/07/29
DAY = "2024-07-29"
ccm.sig_old <- readRDS(paste('ccm_sig', da.name,'nin', nin, DAY, '.obj', sep = '_'))
ccm.rho_old <- readRDS(paste('ccm_rho', da.name,'nin', nin, DAY, 'demo.obj', sep = '_'))
ccm.sig
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 1 1 1
[3,] 1 1 1 1
[4,] 1 1 1 1
ccm.rho
[,1] [,2] [,3] [,4]
[1,] 0.9968866 0.4473342 0.9284430 0.5018034
[2,] 0.3379242 0.9957334 0.5413753 0.8980526
[3,] 0.9220910 0.4914316 0.9979308 0.4655782
[4,] 0.4084419 0.9041382 0.3693686 0.9952968
Based on the causal variables detected in CCM,
generates the time lag series (with maximum lag = 3),
conduct multivariate simplex projection and evaluate the forecast skill, and
selected the best kn(=100) SSR.
Perform multiview embedding analysis for each node
Warning: It is time consuming for running multiview embedding for each nodes
esele_lag <- esim.lag.demo(ds, ccm.rho, ccm.sig, Ed, kmax = 10000, kn = 100, max_lag = 3, Emax = Emax)
sample 1 complete, 37.721 sec
sample 2 complete, 38.177 sec
sample 3 complete, 38.148 sec
sample 4 complete, 37.716 sec
# To avoid overwrite the original files, we save them with different names, 'XXX_NEW'.
DAY = Sys.Date()
saveRDS(esele_lag, paste('eseleLag', da.name, 'nin', nin, DAY, 'demo.obj', sep='_'))
The list of the selected best multivariate SSRs for each node (species, sp) with
the correlation coefficient (rho)
the list of variables (X_1, X_2, …) used in each SSR, noting that the number of variables is identical to the optimal embedding dimension obtained in [2.1], and
the list of time lags (vlag_1, vlag_2, …) used in the corresponding variables (X_1, X_2, …).
as.data.frame(esele_lag)
dmatrix.mv <- mvdist.demo(ds, ds.all, esele_lag)
dmatrix.train.mvx <- dmatrix.mv[['dmatrix.train.mvx']]
dmatrix.test.mvx <- dmatrix.mv[['dmatrix.test.mvx']]
Leave-one-out cross-validation for finding the optimal parameters for MDR S-map analysis
–Warning: The cross-validation is the most time-consuming step in MDR S-map requiring massive computations and . Thus, we recommend dividing job into smaller parts (sub.da>1) or used parallel computation (parall=T, ncore>=1)
do.MDR.CV <- T
### The parameter cv.unit determines the precision of selected parameters and strongly influences computation time.
### In our cases, we used cv.unit=0.025 to obtain more precise estimations
### This parameter may be adjusted to 0.05 or even 0.1, depending on how sensitive the results to parameter precision.
cv.unit <- 0.025
alpha.so <- seq(0, 1, cv.unit); # Sequence of alpha
sub.da <- 1 # Divide the computation job into five parts
afsp <- eqsplit(1:length(alpha.so), sub.da) # Divide the parameter space based on alpha parameter
alf <- 1 # Run CV in the first parameter subset
# Cross-validation of MDR analysis
if(do.MDR.CV){
alpha.s <- alpha.so[afsp[alf, 1]:afsp[alf, 2]] # Subset parameter space
cv.ind <- cv.MDR.demo(ds, ds_tp1, dmatrix.list = dmatrix.train.mvx,
parall = T, ncore = 24, keep_intra = T,alpha.seq = alpha.s)
# To avoid overwrite the original files, we save them with different names, 'XXX_NEW'.
saveRDS(cv.ind, paste(da.name, 'nin', nin, 'cvunit', cv.unit, 'alph', alpha.s[1]*100, DAY, 'cvout_Nmvx_Rallx.obj', sep = '_'))
}
Select the optimal parameter set for regularized regression model with the minimal MSE
paracv.demo <- secv.demo(cv.ind)
paracv.demo
Fitting MDR S-map based on the parameters selected by CV
Note that the linear regression model in MDR S-map includes all of variables as X even if some of them were not selected as the causal variables at the step of using CCM. This is reasonable because S-map coefficients represent the effects of variables at short timescale while CCM evaluates the causality in the whole period of the time series.
#Setting
do.MDR <- F
cv.unit <- 0.025
ptype <- 'aenet' #enet:elastic-net or msaenet: adaptive elastic-net
#Fitting the MDR S-map
smap.2H2P <- MDRsmap.demo(paracv = paracv.demo, ptype = ptype, keep_intra = T, out.sample = T,
ds,ds_tp1,ds.test,dst_tp1,
dmatrix.list = dmatrix.train.mvx,
dmatrix.test.list = dmatrix.test.mvx)
In case for saving the results
DAY = Sys.Date()
nr.out <- smap.2H2P[['nr.out']];
saveRDS(nr.out, paste(da.name,'_nin', nin, 'cvunit', cv.unit, ptype, DAY, 'nrout_Nmvx_Rallx_demo_NEW.obj',sep = '_'))
# Save interaction Jacobian matrices at all time points
saveRDS(smap.2H2P[['jcof']], paste(da.name,'nin', nin, 'cvunit', cv.unit, ptype, DAY, '_jcof_Nmvx_Rallx_demo_NEW.obj', sep='_'))
S-map coefficients
smap.2H2P[['jcof']]
Loading data
smap.2H2P <- list()
smap.2H2P[['jcof']] <- readRDS("4p_model_nin_198_cvunit_0.025_aenet_2024-07-31__jcof_Nmvx_Rallx_demo_NEW.obj")
smap.2H2P[['jcof']]
final_ID <- 200
delta_11 <- 1
MDR_Smap_dH1dH1 <- subset(smap.2H2P[['jcof']], variable == j_Host1)[1:199, 5] - delta_11
MDR_Smap_dH1dH2 <- subset(smap.2H2P[['jcof']], variable == j_Host1)[1:199, 6]
MDR_Smap_dH1dP1 <- subset(smap.2H2P[['jcof']], variable == j_Host1)[1:199, 7]
MDR_Smap_dH1dP2 <- subset(smap.2H2P[['jcof']], variable == j_Host1)[1:199, 8]
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11, type = "l", col = 1, ylim = c(-0.9, 1.0), xlab = "time", ylab = "Interation coefficients")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_2[-final_ID], type = "l", col = "red", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID], type = "l", col = "blue", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_4[-final_ID], type = "l", col = "green", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], MDR_Smap_dH1dH1, type = "l", col = "brown", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
par(new = T)
plot(HP_result$t[-final_ID], MDR_Smap_dH1dP1, type = "l", col = "cyan", ylim = c(-0.9, 1.0), xlab = "", ylab = "")
legend("topright",
legend = c("S-map from H1 to H1 - delta11", "S-map from H2 to H1", "S-map from P1 to H1", "S-map from P2 to H1","MDR S-map from H1 to H1 - delta11", "MDR S-map from P1 to H1"),
col = c("black", "red", "blue", "green", "brown", "cyan"),
lty = c(1, 1, 1, 1, 1, 1),
cex = 0.8)
\[ \frac{\partial F_{H,i}}{\partial H_{i,k}} = -1 + (1 - a_i H_i) e ^ {-a_i H_i - c_{1i} P_1 - c_{2i} P_2} (r_i + e ^{-m_{H,i}}) \] \[ \frac{\partial F_{H,i}}{\partial H_{i',k}} = 0 \] \[ \frac{\partial F_{H,i}}{\partial P_{j,k}} = c_{ji} e ^ {-a_i H_i - c_{1i} P_1 - c_{2i} P_2} (-r_i - e ^ {-m_{H,i}}) H_k \]
Functions for partial derivatives for H1:
dF_H1dH1 <- function(data) {
coeff <- -1 + (1 - a1 * data$H1) * exp(- a1 * data$H1 - c11 * data$P1 - c21 * data$P2) * (r1 + exp(- mH1))
return(coeff)
}
dF_H1dP1 <- function(data) {
coeff <- c11 * exp(- a1 * data$H1 - c11 * data$P1 - c21 * data$P2) * (- r1 - exp(- mH1)) * data$H1
return(coeff)
}
dF_H1dP2 <- function(data) {
coeff <- c21 * exp(- a1 * data$H1 - c11 * data$P1 - c21 * data$P2) * (- r1 - exp(- mH1)) * data$H1
return(coeff)
}
Dataframe for ggplot
#For lineplots
H2P2_H1toH1_ggplot <- data.frame(
time = HP_result$t[-final_ID],
IntS = HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11,
type = rep("01std_smapH1toH1", 199)
)
H2P2_H1toH1_ggplot <- rbind.data.frame(
H2P2_H1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = MDR_Smap_dH1dH1,
type = rep("02mdr_smapH1toH1", 199)
)
)
H2P2_H1toH1_ggplot <- rbind.data.frame(
H2P2_H1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = dF_H1dH1(HP_result)[-final_ID],
type = rep("03IS_dH1dH1", 199)
)
)
H2P2_P1toH1_ggplot <- data.frame(
time = HP_result$t[-final_ID],
IntS = HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID],
type = rep("01std_smapP1toH1", 199)
)
H2P2_P1toH1_ggplot <- rbind.data.frame(
H2P2_P1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = MDR_Smap_dH1dP1,
type = rep("02mdr_smapP1toH1", 199)
)
)
H2P2_P1toH1_ggplot <- rbind.data.frame(
H2P2_P1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1],
type = rep("03IS_dH1dP1", 199)
)
)
#For error plots
H2P2_error_H1toH1_ggplot <- data.frame(
time = HP_result$t[-final_ID],
IntS = HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11 - dF_H1dH1(HP_result)[-final_ID],
type = rep("01std_smapH1toH1", 199)
)
H2P2_error_H1toH1_ggplot <- rbind.data.frame(
H2P2_error_H1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = MDR_Smap_dH1dH1 - dF_H1dH1(HP_result)[-final_ID],
type = rep("02mdr_smapH1toH1", 199)
)
)
H2P2_error_P1toH1_ggplot <- data.frame(
time = HP_result$t[-final_ID],
IntS = HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID] - dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1],
type = rep("01std_smapP1toH1", 199)
)
H2P2_error_P1toH1_ggplot <- rbind.data.frame(
H2P2_error_P1toH1_ggplot, data.frame(
time = HP_result$t[-final_ID],
IntS = MDR_Smap_dH1dP1 - dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1],
type = rep("02mdr_smapP1toH1", 199)
)
)
fig2a_ggplot <- ggplot(data = H2P2_H1toH1_ggplot) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00")) +
scale_shape_manual(values = c(15, 16, 17)) +
 geom_line(aes(x = time, y = IntS, group = type, color = type)) +
geom_point(aes(x = time, y = IntS, color = type, shape = type), size = 2) +
xlim(1800, 1999) + theme_bw() + theme(legend.position = "bottom") +
labs(title = "Figure 2a: Time evolution of H1 -> H1", x = "time", y = "H1 -> H1")
fig2b_ggplot <- ggplot(H2P2_error_H1toH1_ggplot, aes(type, abs(IntS))) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00")) +
scale_shape_manual(values = c(15, 16, 17)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(shape = type, color = type, alpha = 0.1), width = 0.25, size = 3) +
theme_bw() + theme(legend.position = "bottom") +
labs (title = "Figure 2b: Comparison to parametric IS: H1 -> H1")
plot(fig2a_ggplot)
ggsave("fig2a.pdf", width = 6, height = 4)
plot(fig2b_ggplot)
ggsave("fig2b.pdf", width = 4, height = 4)
fig2c_ggplot <- ggplot(data = H2P2_P1toH1_ggplot) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00")) +
scale_shape_manual(values = c(15, 16, 17)) +
 geom_line(aes(x = time, y = IntS, group = type, color = type)) +
geom_point(aes(x = time, y = IntS, color = type, shape = type), size = 2) +
xlim(1800, 1999) + theme_bw() + theme(legend.position = "bottom") +
labs(title = "Figure 2c: Time evolution of P1 -> H1", x = "time", y = "P1 -> H1")
fig2d_ggplot <- ggplot(H2P2_error_P1toH1_ggplot, aes(type, abs(IntS))) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00")) +
scale_shape_manual(values = c(15, 16, 17)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(shape = type, color = type, alpha = 0.1), width = 0.25, size = 3) +
theme_bw() + theme(legend.position = "bottom") +
labs (title = "Figure 2d: Comparison to parametric IS: P1 -> H1")
plot(fig2c_ggplot)
ggsave("fig2c.pdf", width = 6, height = 4)
plot(fig2d_ggplot)
ggsave("fig2d.pdf", width = 4, height = 4)
#for diagonal element (H1 -> H1)
plot(dF_H1dH1(HP_result)[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11,
xlab = "theoretical coefficient (H1 -> H1)",
ylab = "standard S-map coefficient (H1 -> H1)",
)
abline(coef = c(0,1))
cor(dF_H1dH1(HP_result)[-final_ID], HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11)
[1] 0.9781177
#Y should be the theoretical coefficient
x <- HP_m_smap_res2_s$smap_coefficients[[1]]$c_1[-final_ID] - delta_11
summary(lm(dF_H1dH1(HP_result)[-final_ID] ~ x))
Call:
lm(formula = dF_H1dH1(HP_result)[-final_ID] ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.082187 -0.024128 0.001541 0.023675 0.097218
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.001131 0.002422 -0.467 0.641
x 1.047491 0.015874 65.986 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03386 on 197 degrees of freedom
Multiple R-squared: 0.9567, Adjusted R-squared: 0.9565
F-statistic: 4354 on 1 and 197 DF, p-value: < 2.2e-16
plot(dF_H1dH1(HP_result)[-final_ID], MDR_Smap_dH1dH1,
xlab = "theoretical coefficient (H1 -> H1)",
ylab = "MDR S-map coefficient (H1 -> H1)"
)
abline(coef = c(0,1))
cor(dF_H1dH1(HP_result)[-final_ID], MDR_Smap_dH1dH1)
[1] 0.6705204
summary(lm(dF_H1dH1(HP_result)[-final_ID] ~ MDR_Smap_dH1dH1))
Call:
lm(formula = dF_H1dH1(HP_result)[-final_ID] ~ MDR_Smap_dH1dH1)
Residuals:
Min 1Q Median 3Q Max
-0.32306 -0.07849 -0.01325 0.07168 0.33871
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.040208 0.009882 4.069 6.83e-05 ***
MDR_Smap_dH1dH1 1.832602 0.144465 12.685 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1208 on 197 degrees of freedom
Multiple R-squared: 0.4496, Adjusted R-squared: 0.4468
F-statistic: 160.9 on 1 and 197 DF, p-value: < 2.2e-16
#for off-diagonal element (P1 -> H1)
plot(dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1], HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID],
xlab = "theoreticalcoefficient (P1 -> H1)",
ylab = "standard S-map coefficient (P1 -> H1)"
)
abline(coef = c(0,1))
cor(dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1], HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID])
[1] 0.9619422
y <- dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1]
summary(lm(y ~ HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID]))
Call:
lm(formula = y ~ HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID])
Residuals:
Min 1Q Median 3Q Max
-0.116567 -0.021643 0.001839 0.021454 0.109281
Coefficients:
Estimate Std. Error t value
(Intercept) -0.007259 0.005382 -1.349
HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID] 1.014493 0.020532 49.410
Pr(>|t|)
(Intercept) 0.179
HP_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID] <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03924 on 197 degrees of freedom
Multiple R-squared: 0.9253, Adjusted R-squared: 0.925
F-statistic: 2441 on 1 and 197 DF, p-value: < 2.2e-16
plot(dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1], MDR_Smap_dH1dP1,
xlab = "theoreticalcoefficient (P1 -> H1)",
ylab = "MDR S-map coefficient (P1 -> H1)"
)
abline(coef = c(0,1))
cor(MDR_Smap_dH1dP1, dF_H1dP1(HP_result)[-final_ID] * HP_result.sd[3] / HP_result.sd[1])
[1] 0.9709203
summary(lm(y ~ MDR_Smap_dH1dP1))
Call:
lm(formula = y ~ MDR_Smap_dH1dP1)
Residuals:
Min 1Q Median 3Q Max
-0.122178 -0.017208 0.004087 0.022169 0.071996
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.049246 0.005555 8.864 4.6e-16 ***
MDR_Smap_dH1dP1 1.368911 0.024048 56.923 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03438 on 197 degrees of freedom
Multiple R-squared: 0.9427, Adjusted R-squared: 0.9424
F-statistic: 3240 on 1 and 197 DF, p-value: < 2.2e-16
\[ \frac{dP_i}{dt} = v_i \lambda_i \frac{P_i C_i}{C_i + C_i^*} - v_i P_i, \quad (i=1,2) \tag{1} \]
\[ \frac{dC_i}{dt} = \mu_i \kappa_i \frac{C_i R}{R + R^*} - v_i \lambda_i \frac{P_i C_i}{C_i + C_i^*} - \mu_i C_i, \quad (i=1,2) \tag{2} \]
\[ \frac{dR}{dt} = R \left( 1 - \frac{R}{k} \right) - \sum_{i=1,2} \mu_i \kappa_i \frac{C_iR}{R + R^*}.\tag{3} \]
v1 <- 0.1
v2 <- 0.07
lambd1 <- 3.2
lambd2 <- 2.9
Cstar1 <- 0.5
Cstar2 <- 0.5
myu1 <- 0.15
myu2 <- 0.15
kp1 <- 2.5
kp2 <- 2.0
Rstar <- 0.3
k <- 1.2
j_P1 <- 1
j_P2 <- 2
j_C1 <- 3
j_C2 <- 4
j_R <- 5
dP1dt <- function(in_vec, t)
{
growth <- v1 * lambd1 * in_vec[j_P1] * in_vec[j_C1] / (in_vec[j_C1] + Cstar1)
mortality <- v1 * in_vec[j_P1]
return(growth - mortality)
}
dP2dt <- function(in_vec, t)
{
growth <- v2 * lambd2 * in_vec[j_P2] * in_vec[j_C2] / (in_vec[j_C2] + Cstar2)
mortality <- v2 * in_vec[j_P2]
return(growth - mortality)
}
dC1dt <- function(in_vec, t)
{
growth <- myu1 * kp1 * in_vec[j_C1] * in_vec[j_R] / (in_vec[j_R] + Rstar)
mortality <- v1 * lambd1 * in_vec[j_P1] * in_vec[j_C1] / (in_vec[j_C1] + Cstar1) + myu1 * in_vec[j_C1]
return(growth - mortality)
}
dC2dt <- function(in_vec, t)
{
growth <- myu2 * kp2 * in_vec[j_C2] * in_vec[j_R] / (in_vec[j_R] + Rstar)
mortality <- v2 * lambd2 * in_vec[j_P2] * in_vec[j_C2] / (in_vec[j_C2] + Cstar2) + myu2 * in_vec[j_C2]
return(growth - mortality)
}
dRdt <- function(in_vec, t)
{
growth <- in_vec[j_R] * (1.0 - in_vec[j_R] / k)
mortality <- myu1 * kp1 * in_vec[j_C1] * in_vec[j_R] / (in_vec[j_R] + Rstar) + myu2 * kp2 * in_vec[j_C2] * in_vec[j_R] / (in_vec[j_R] + Rstar)
return(growth - mortality)
}
#function to calculate all coefficients, five dimensional
diff_5sp <- function(in_vec, ref_data, t, h_interval, dim){
temp_vec <- numeric(dim)
temp_vec[j_P1] <- dP1dt(in_vec, t)
temp_vec[j_P2] <- dP2dt(in_vec, t)
temp_vec[j_C1] <- dC1dt(in_vec, t)
temp_vec[j_C2] <- dC2dt(in_vec, t)
temp_vec[j_R] <- dRdt(in_vec, t)
temp_vec
}
dim_model1 <- 5
nv0 <- numeric(dim_model1)
nv <- numeric(dim_model1)
nv0[j_R] <- 1.0
nv0[j_C1] <- 0.5
nv0[j_C2] <- 0.8
nv0[j_P1] <- 0.7
nv0[j_P2] <- 0.8
deltat <- 0.01
t <- 0.0 # initial condition (initial time, 0)
end_time <- 2000; tau <- 5; transient_period <- 1000
write_index <- 1
#Initial condition set
nv <- nv0
#write initial condition
cat(t, nv, "\n")
0 0.7 0.8 0.5 0.8 1
#For the transient dynamics (RUN BURN)
for(i in 1:as.integer(transient_period / deltat)){
nv <- rk4(in_vec = nv, time = t, h_interval = deltat, dim = dim_model1, diff_vec = diff_5sp)
t <- t + deltat # update time
}
# record the initial condition after BURNING phase
rk4_result <- data.frame(t = t, P1 = nv[j_P1], P2 = nv[j_P2], C1 = nv[j_C1], C2 = nv[j_C2], R = nv[j_R])
#After transient period
for(i in 1:as.integer((end_time - transient_period) / deltat)){
nv <- rk4(in_vec = nv, time = t, h_interval = deltat, dim = dim_model1, diff_vec = diff_5sp)
t <- t + deltat # update time
###record the result every "tau only
if(write_index < as.integer(tau / deltat)){
write_index <- write_index + 1
}
else {
#cat(time, nv_ee1,"\n")
rk4_result <- rbind.data.frame(rk4_result, c(t, nv[j_P1], nv[j_P2], nv[j_C1], nv[j_C2], nv[j_R]))
write_index <- 1
}
}
head(rk4_result)
saveRDS(rk4_result, "rk4_result_5sp_model.obj")
plot(rk4_result$t, rk4_result$C1, col = "red", type = "l", xlab = "time", ylab = "abundance", ylim = c(0, 2.0))
par(new = T)
plot(rk4_result$t, rk4_result$C2, col = "blue", type = "l", xlab = "", ylab = "", ylim = c(0, 2.0))
par(new = T)
plot(rk4_result$t, rk4_result$R, col = "green", type = "l", xlab = "", ylab = "", ylim = c(0, 2.0))
The partial derivative of dC1/dt
with respect to C1: \[ \frac{\partial}{\partial C_{1}} \left( \mu_{1} K_{1} \frac{C_{1} R}{R + R^{*}} - v_{1} \lambda_{1} \frac{P_{1} C_{1}}{C_{1} + C_{1}^{*}} - \mu_{1} C_{1} \right) = K_{1} \mu_{1} \frac{R}{R + R^{*}} - v_{1} \lambda_{1} \frac{P_{1} C_{1}*}{(C_{1} + C_{1}^{*})^2} - \mu_{1} \]
with respect to R: \[ \frac{\partial}{\partial R} \left( \mu_{1} \kappa_{1} \frac{C_{1} R}{R + R^{*}} \right) = \kappa_{1} \mu_{1} \frac{C_{1} R^*}{(R + R^{*})^2} \]
Function definition
dC1dC1 <- function(data) {
coeff <- myu1 * kp1 * data$R / (data$R + Rstar) - v1 * lambd1 * data$P1 * Cstar1/ (data$C1 + Cstar1)^2 - myu1
return(coeff)
}
dC1dR <- function(data) {
coeff <- myu1 * kp1 * data$C1 *Rstar / (data$R + Rstar)^2
return(coeff)
}
Based on https://ushio-ecology-blog.blogspot.com/2019/12/20191225blogger0007.html
rk4_result.mean <- apply(rk4_result[, -1], 2, mean, na.rm = T) # mean abundance
rk4_result.sd <- apply(rk4_result[, -1], 2, sd, na.rm = T) # SD of abundance
mean_mat <- matrix(rk4_result.mean, nrow = nrow(rk4_result) , ncol = length(rk4_result.mean), byrow = TRUE)
sd_mat <- matrix(rk4_result.sd, nrow = nrow(rk4_result) , ncol = length(rk4_result.mean), byrow = TRUE)
rk4_result_s <- (rk4_result[, -1] - mean_mat) /sd_mat # Standardized dataset
head(rk4_result_s)
# data block for smap
rk4_smap_block_s <- data.frame(P1 = rk4_result_s$P1, P2 = rk4_result_s$P2, C1 = rk4_result_s$C1, C2 = rk4_result_s$C2, R = rk4_result_s$R)
# delete a few column
#rk4_smap_block_s <- rk4_smap_block_s[, -j_P2] #delete P2 from the embedding
# conducting smap for obtaining the optimal theta
rk4_m_smap_res_s <- block_lnlp(rk4_smap_block_s, method = "s-map",
target_column = j_C1,
theta = c(0, 1e-04, 3e-04, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8),
silent = TRUE)
Plotting the effect of theta
plot(rk4_m_smap_res_s$theta, rk4_m_smap_res_s$rmse, type = "b")
(rk4_opt_theta_s <- rk4_m_smap_res_s$theta[which.min(rk4_m_smap_res_s$rmse)])
[1] 8
Conducting standard smap with the optimal theta
rk4_m_smap_res2_s <- block_lnlp(rk4_smap_block_s, method = "s-map",
target_column = j_C1,
theta = rk4_opt_theta_s,
silent = TRUE,
save_smap_coefficients = TRUE)
head(as.data.frame(rk4_m_smap_res2_s$smap_coefficients[[1]]))
final_ID <- 201
plot(rk4_result$t[-final_ID], rk4_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID] - 1.0, type = "l", col = 1, ylim = c(-1, 1), xlab = "time", ylab = "Smap coefficient: C1 -> C1")
par(new = T)
plot(rk4_result$t[-final_ID], tau*dC1dC1(rk4_result)[-final_ID], type = "l", col = "blue", ylim = c(-1, 1), xlab = "", ylab = "")
\[ \frac{\partial F_{P_1}}{\partial P_1} = v_1 \lambda_1 \frac{C_1}{C_1 + C_1^*} - v_1, \frac{\partial F_{P_1}}{\partial P_2} = 0, \frac{\partial F_{P_1}}{\partial C_1} = v_1 \lambda_1 \frac{P_1 C_1^*}{(C_1 + C_1^*)^2}, \frac{\partial F_{P_1}}{\partial C_2} = \frac{\partial F_{P_1}}{\partial R} = 0, \]
\[ \frac{\partial F_{P_2}}{\partial P_1} = 0, \frac{\partial F_{P_2}}{\partial P_2} = v_2 \lambda_2 \frac{C_2}{C_2 + C_2^*} - v_2, \frac{\partial F_{P_2}}{\partial C_1} = 0, \frac{\partial F_{P_2}}{\partial C_2} = v_2 \lambda_2 \frac{P_2 C_2^*}{(C_2 + C_2^*)^2}, \frac{\partial F_{P_2}}{\partial R} = 0, \]
\[ \frac{\partial F_{C_1}}{\partial P_1} = -v_1 \lambda_1 \frac{C_1}{C_1 + C_1^*}, \frac{\partial F_{C_1}}{\partial P_2} = 0, \frac{\partial F_{C_1}}{\partial C_1} = \mu_1 \kappa_1 \frac{R}{R + R^*} - v_1 \lambda_1 \frac{P_1 C_1^*}{(C_1 + C_1^*)^2} - \mu_1, \frac{\partial F_{C_1}}{\partial C_2} = 0, \frac{\partial F_{C_1}}{\partial R} = \mu_1 \kappa_1 \frac{C_1 R^*}{(R + R^*)^2}, \]
\[ \frac{\partial F_{C_2}}{\partial P_1} = 0, \frac{\partial F_{C_2}}{\partial P_2} = -v_2 \lambda_2 \frac{C_2}{C_2 + C_2^*}, \frac{\partial F_{C_2}}{\partial C_1} = 0, \frac{\partial F_{C_2}}{\partial C_2} = \mu_2 \kappa_2 \frac{R}{R + R^*} - v_2 \lambda_2 \frac{P_2 C_2^*}{(C_2 + C_2^*)^2} - \mu_2, \frac{\partial F_{C_2}}{\partial R} = \mu_2 \kappa_2 \frac{C_2 R^*}{(R + R^*)^2}, \]
\[ \frac{\partial F_{R}}{\partial P_1} = \frac{\partial F_{R}}{\partial P_2} = 0, \frac{\partial F_{R}}{\partial C_1} = - \mu_1 \kappa_1 \frac{R}{R + R^*}, \frac{\partial F_{R}}{\partial C_2} = - \mu_2 \kappa_2 \frac{R}{R + R^*}, \frac{\partial F_{R}}{\partial R} = \left( 1 - \frac{2 R}{k} \right)- \sum_{i=1,2} \mu_i \kappa_i \frac{C_i R^*}{(R + R^*)^2}. \]
Partial derivatives
rF_P1rX <- function(in_vec, t) {
nv_p <- numeric(length(in_vec))
nv_p[j_P1] <- v1 * lambd1 * in_vec[j_C1] / (in_vec[j_C1] + Cstar1) - v1
nv_p[j_P2] <- 0
nv_p[j_C1] <- v1 * lambd1 * in_vec[j_P1] * Cstar1 / (in_vec[j_C1] + Cstar1)^2
nv_p[j_C2] <- 0
nv_p[j_R] <- 0
return(nv_p)
}
rF_P2rX <- function(in_vec, t) {
nv_p <- numeric(length(in_vec))
nv_p[j_P2] <- v2 * lambd2 * in_vec[j_C2] / (in_vec[j_C2] + Cstar2) - v2
nv_p[j_P1] <- 0
nv_p[j_C2] <- v2 * lambd2 * in_vec[j_P2] * Cstar2 / (in_vec[j_C2] + Cstar2)^2
nv_p[j_C1] <- 0
nv_p[j_R] <- 0
return(nv_p)
}
rF_C1rX <- function(in_vec, t) {
nv_p <- numeric(length(in_vec))
nv_p[j_P1] <- -v1 * lambd1 * in_vec[j_C1] / (in_vec[j_C1] + Cstar1)
nv_p[j_P2] <- 0
nv_p[j_C1] <- myu1 * kp1 * in_vec[j_R] / (in_vec[j_R] + Rstar) - v1 * lambd1 * in_vec[j_P1] * Cstar1/ (in_vec[j_C1] + Cstar1)^2 - myu1
nv_p[j_C2] <- 0
nv_p[j_R] <- myu1 * kp1 * in_vec[j_C1] * Rstar / (in_vec[j_R] + Rstar)^2
return(nv_p)
}
rF_C2rX <- function(in_vec, t) {
nv_p <- numeric(length(in_vec))
nv_p[j_P2] <- -v2 * lambd2 * in_vec[j_C2] / (in_vec[j_C2] + Cstar2)
nv_p[j_P1] <- 0
nv_p[j_C2] <- myu2 * kp2 * in_vec[j_R] / (in_vec[j_R] + Rstar) - v2 * lambd2 * in_vec[j_P2] * Cstar2/ (in_vec[j_C2] + Cstar2)^2 - myu2
nv_p[j_C1] <- 0
nv_p[j_R] <- myu2 * kp2 * in_vec[j_C2] * Rstar / (in_vec[j_R] + Rstar)^2
return(nv_p)
}
rF_RrX <- function(in_vec, t) {
nv_p <- numeric(length(in_vec))
nv_p[j_P1] <- 0
nv_p[j_P1] <- 0
nv_p[j_C1] <- -myu1 * kp1 * in_vec[j_R] / (in_vec[j_R] + Rstar)
nv_p[j_C2] <- -myu2 * kp2 * in_vec[j_R] / (in_vec[j_R] + Rstar)
nv_p[j_R] <- (1 - 2 * in_vec[j_R] / k) - (myu1 * kp1 * in_vec[j_C1] + myu2 * kp2 * in_vec[j_C2]) * Rstar / (in_vec[j_R] + Rstar)^2
return(nv_p)
}
Linearlized perturbed ODE system equations \[ \frac{d\boldsymbol{\delta}(t)}{dt} = \mathbf{J}_{t} \boldsymbol{\delta}(t), \quad \boldsymbol{\delta}(0) = \Delta \boldsymbol{x}_{j}. \tag{7a} \]
\[\begin{align} IC_{i,j,adj}(k=k_1) &= \lim_{\Delta x \to 0} \frac{\left[\phi(\boldsymbol{x}(k_1) + \Delta x_j, \tau) - \phi(\boldsymbol{x}(k_1), \tau)\right]_i}{\Delta x} - \delta_{ij}\\ &= \lim_{\Delta x \to 0} \frac{\left[\boldsymbol{\delta}_n (\tau)\right]_i}{\Delta x} - \delta_{ij} \\ &\approx \frac{\left[\boldsymbol{\delta}_n (\tau)\right]_i}{\Delta x} - \delta_{ij} \tag{7c} \end{align}\]
Linearlized ODE system
dP1Ldt <- function(in_vec, ref_data, t)
{
ref_vec <- as.numeric(ref_data[t, -1])
return(as.numeric(rF_P1rX(ref_vec) %*% in_vec))
}
dP2Ldt <- function(in_vec, ref_data, t)
{
ref_vec <- as.numeric(ref_data[t, -1])
return(as.numeric(rF_P2rX(ref_vec) %*% in_vec))
}
dC1Ldt <- function(in_vec, ref_data, t)
{
ref_vec <- as.numeric(ref_data[t, -1])
return(as.numeric(rF_C1rX(ref_vec) %*% in_vec))
}
dC2Ldt <- function(in_vec, ref_data, t)
{
ref_vec <- as.numeric(ref_data[t, -1])
return(as.numeric(rF_C2rX(ref_vec) %*% in_vec))
}
dRLdt <- function(in_vec, ref_data, t)
{
ref_vec <- as.numeric(ref_data[t, -1])
return(as.numeric(rF_RrX(ref_vec) %*% in_vec))
}
#function to calculate all coefficients, five dimensional
diff_L_5sp <- function(in_vec, ref_data, t, h_interval, dim){
temp_vec <- numeric(dim)
t_index <- as.integer((t / h_interval) * 2) #to convert continuous time to the number of row
temp_vec[j_P1] <- dP1Ldt(in_vec, ref_data, t_index)
temp_vec[j_P2] <- dP2Ldt(in_vec, ref_data, t_index)
temp_vec[j_C1] <- dC1Ldt(in_vec, ref_data, t_index)
temp_vec[j_C2] <- dC2Ldt(in_vec, ref_data, t_index)
temp_vec[j_R] <- dRLdt(in_vec, ref_data, t_index)
temp_vec
}
Since rk4() includes the evaluation of derivatives at t + deltat/2, we need to have numerical solutions with the high resolution with deltat/2
t <- 0.0 # initial condition (initial time, 0)
end_time <- 2000; tau <- 5; transient_period <- 1000
write_index <- 1
deltat_fine <- deltat * 0.5
#Initial condition set
nv_rk4_fine <- nv0
# write initial condition
cat(t, nv_rk4_fine, "\n")
0 0.7 0.8 0.5 0.8 1
#For the transient dynamics (RUN BURN)
for(i in 1:as.integer(transient_period / deltat_fine)){
nv_rk4_fine <- rk4(in_vec = nv_rk4_fine, time = t, h_interval = deltat_fine, dim = dim_model1, diff_vec = diff_5sp)
t <- t + deltat_fine # update time
}
# record the initial condition after BURNING period
rk4_result_fine <- data.frame(t = t, P1 = nv_rk4_fine[j_P1], P2 = nv_rk4_fine[j_P2], C1 = nv_rk4_fine[j_C1], C2 = nv_rk4_fine[j_C2], R = nv_rk4_fine[j_R])
#After transient period
for(i in 1:as.integer((end_time - transient_period) / deltat_fine)){
nv_rk4_fine <- rk4(in_vec = nv_rk4_fine, time = t, h_interval = deltat_fine, dim = dim_model1, diff_vec = diff_5sp)
t <- t + deltat_fine # update time
rk4_result_fine <- rbind.data.frame(rk4_result_fine, c(t, nv_rk4_fine[j_P1], nv_rk4_fine[j_P2], nv_rk4_fine[j_C1], nv_rk4_fine[j_C2], nv_rk4_fine[j_R]))
}
head(rk4_result_fine)
#saveRDS(rk4_result_fine, paste("default", "_rk4_result_fine.obj", sep=""))
#saveRDS(rk4_result_fine, paste(Sys.Date(), "_rk4_result_fine.obj", sep=""))
Comparing the solution
plot(rk4_result$t, rk4_result$C1, col = "red", type = "l", xlab = "time", ylab = "abundance", ylim = c(0, 2.0))
par(new = T)
plot(rk4_result_fine$t, rk4_result_fine$C1, col = "black", type = "l", xlab = "", ylab = "", ylim = c(0, 2.0), lty = "dashed")
#loading the default
rk4_result_fine <- readRDS("default_rk4_result_fine.obj")
delta_x = 0.01 #size of a small perturbation
t <- deltat / 2 # initial condition (initial time, t_abs = transient_period)
t_abs <- t + transient_period - deltat / 2
end_time <- 1.0; tau <- 5;
write_index <- 1
#Initial condition
nv_linear_p <- c(0, 0, delta_x, 0, 0)
# record the initial condition
result_linear_p <- data.frame(t = t_abs, P1 = nv_linear_p[j_P1], P2 = nv_linear_p[j_P2], C1 = nv_linear_p[j_C1], C2 = nv_linear_p[j_C2], R = nv_linear_p[j_R])
for(i in 1:as.integer(end_time / deltat)){
nv_linear_p <- rk4(in_vec = nv_linear_p, ref_data = rk4_result_fine, time = t, h_interval = deltat, dim = dim_model1, diff_vec = diff_L_5sp)
t <- t + deltat # update time (starting from the middle)
t_abs <- t_abs + deltat
result_linear_p <- rbind.data.frame(result_linear_p, c(t_abs, nv_linear_p[j_P1], nv_linear_p[j_P2], nv_linear_p[j_C1], nv_linear_p[j_C2], nv_linear_p[j_R]))
}
result_linear_p
NA
The direct evaluation by the nonlinear ODE
t <- 1000 - transient_period + deltat / 2 # initial condition (initial time, t_abs = start_t)
t_index_start <- as.integer((t / deltat) * 2)
t_index_end <- t_index_start + as.integer((1.0 / deltat) * 2)
#Initial condition
nv_nonlinear_p <- as.numeric(rk4_result_fine[t_index_start, -1]) + c(0, 0, delta_x, 0, 0)
# record the initial condition
result_nonlinear_p <- data.frame(t = t_abs, P1 = nv_nonlinear_p[j_P1], P2 = nv_nonlinear_p[j_P2], C1 = nv_nonlinear_p[j_C1], C2 = nv_nonlinear_p[j_C2], R = nv_nonlinear_p[j_R])
for(i in 1:as.integer(end_time / deltat)){
nv_nonlinear_p <- rk4(in_vec = nv_nonlinear_p, time = t, h_interval = deltat, dim = dim_model1, diff_vec = diff_5sp)
t <- t + deltat # update time (starting from the middle)
}
nv_nonlinear_p - as.numeric(rk4_result_fine[t_index_end, -1]) #difference between perturbed and unperturbed solutions
[1] 4.152205e-05 -4.908053e-13 1.139558e-02 -5.768626e-09 -2.217930e-03
#ref_data: the dataset that includes the fine resolution solution of the target nonlinear ODE systems, of which record started at absolute time = transient_period
#deltax: the size of perturbation at t = h_interval/2
#start_t: usually specified as transient_period
#tau: the interval for the calculation
#dim: dimension of ODE system
#diff_vec: the vector field generated by the linearized ODEs with non-perturbed solution
precise_coeff <- function(ref_data, deltax, start_t, tau, h_interval, dim, diff_vec) {
#matrix to stock results
nv_linear_p <- matrix(0, nrow = dim, ncol = dim)
#Initial condition
nv_linear_p[, j_P1] <- c(deltax, 0, 0, 0, 0)
nv_linear_p[, j_P2] <- c(0, deltax, 0, 0, 0)
nv_linear_p[, j_C1] <- c(0, 0, deltax, 0, 0)
nv_linear_p[, j_C2] <- c(0, 0, 0, deltax, 0)
nv_linear_p[, j_R] <- c(0, 0, 0, 0, deltax)
tzero <- start_t - transient_period + h_interval / 2 # initial condition (initial time, t_abs = start_t)
end_time <- tau
for(j in 1:dim) {
t <- tzero #initialize time for each of different initial perturbations
for(i in 1:as.integer(end_time / deltat)){
nv_linear_p[, j] <- rk4(nv_linear_p[, j], ref_data, t, h_interval, dim, diff_vec)
t <- t + h_interval # update time (starting from the middle)
}
}
return(nv_linear_p/deltax)
}
#This is a function that directly used the nonlinear ODE with the perturbed initial conditions
direct_coeff <- function(ref_data, deltax, start_t, tau, h_interval, dim, diff_vec) {
#matrix to stock results
nv_nonlinear_p <- matrix(0, nrow = dim, ncol = dim)
tzero <- start_t - transient_period + h_interval / 2 # initial condition (initial time, t_abs = start_t)
t_index_start <- as.integer((tzero / h_interval) * 2)
t_index_end <- t_index_start + as.integer((tau / h_interval) * 2)
#Initial condition
nv_nonlinear_p[, j_P1] <- as.numeric(ref_data[t_index_start, -1]) + c(deltax, 0, 0, 0, 0)
nv_nonlinear_p[, j_P2] <- as.numeric(ref_data[t_index_start, -1]) + c(0, deltax, 0, 0, 0)
nv_nonlinear_p[, j_C1] <- as.numeric(ref_data[t_index_start, -1]) + c(0, 0, deltax, 0, 0)
nv_nonlinear_p[, j_C2] <- as.numeric(ref_data[t_index_start, -1]) + c(0, 0, 0, deltax, 0)
nv_nonlinear_p[, j_R] <- as.numeric(ref_data[t_index_start, -1]) + c(0, 0, 0, 0, deltax)
end_time <- tau
#rk4(in_vec = nv_rk4, time = t, h_interval = deltat, dim = dim_model1, diff_vec = diff_5sp)
for(j in 1:dim) {
t <- tzero
for(i in 1:as.integer(end_time / deltat)){
nv_nonlinear_p[, j] <- rk4(nv_nonlinear_p[, j], ref_data = NULL, t, h_interval, dim, diff_vec)
t <- t + h_interval # update time (starting from the middle)
}
}
#Calculate the difference between the perturbed flow and unperturbed flow
for(j in 1:dim) nv_nonlinear_p[, j] <- nv_nonlinear_p[, j] - as.numeric(ref_data[t_index_end, -1])
return(nv_nonlinear_p/deltax)
}
The reference time points are every tau (= 5.0), but the precise_coeff was calculated for many scenarios (5.0, 0.01, 0.5, 1.0, 2.0, 10.0)
precise_coeff_5sp <- list()
inst_precise_coeff_5sp <- list()
start_time <- vector()
tau <- 5
tau1 <- 5
precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau1, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
tau2 <- 0.01
inst_precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau2, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
tau3 <- 0.5
tau0.5_precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau3, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
tau4 <- 1.0
tau1.0_precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau4, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
tau5 <- 2.0
tau2.0_precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau5, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
tau6 <- 10.0
tau10.0_precise_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_linear <- precise_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau6, deltat, dim = 5, diff_L_5sp)
result_linear
}, mc.cores = 16)
precise_coeff_5sp[[5]]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.656280e-01 5.598464e-07 2.476259e-01 -0.0185625139 8.142788e-01
[2,] 2.989538e-10 7.047559e-01 -1.358904e-09 0.0001775428 5.175086e-08
[3,] -3.732702e-01 3.687060e-06 2.609733e-01 -0.0990582164 3.229635e+00
[4,] 5.365938e-06 -6.729101e-05 -1.929607e-05 0.6277554416 5.591555e-04
[5,] 2.992122e-01 3.685980e-05 -7.940109e-01 -0.7522821462 1.692837e+01
inst_precise_coeff_5sp[[5]]
[,1] [,2] [,3] [,4] [,5]
[1,] 1.000623e+00 4.063713e-18 5.005509e-04 -6.380349e-11 1.479335e-06
[2,] 4.173363e-21 9.993005e-01 -1.026859e-17 6.640377e-07 1.905640e-13
[3,] -1.622704e-03 3.231007e-14 9.981620e-01 -3.814216e-07 5.899615e-03
[4,] 2.510897e-14 -2.537710e-07 -4.635541e-11 9.986299e-01 5.735926e-07
[5,] 1.318010e-07 1.649286e-11 -1.621770e-04 -1.297720e-04 1.003887e+00
tau0.5_precise_coeff_5sp[[5]]
[,1] [,2] [,3] [,4] [,5]
[1,] 1.028152e+00 2.811240e-11 2.607848e-02 -8.868742e-06 4.074820e-03
[2,] 2.589855e-14 9.656173e-01 -1.260250e-12 3.102685e-05 4.676182e-10
[3,] -7.628464e-02 3.957674e-09 9.089931e-01 -9.428815e-04 2.912958e-01
[4,] 3.254568e-09 -1.166418e-05 -1.185950e-07 9.343543e-01 2.943598e-05
[5,] 3.857781e-04 4.751814e-08 -9.475474e-03 -7.675132e-03 1.234747e+00
#saveRDS(precise_coeff_5sp, "default_precise_coeff_5sp.obj")
direct_coeff_5sp <- list()
start_time <- vector()
transient_period <- 1000
direct_coeff_5sp <- mclapply(1:200, function(i) {
start_time[i] <- transient_period + (i - 1)*tau
result_nonlinear <- direct_coeff(rk4_result_fine, deltax = 0.0001, start_t = start_time[i], tau, deltat, dim = 5, diff_5sp)
result_nonlinear
}, mc.cores = 16)
direct_coeff_5sp[[5]]
[,1] [,2] [,3] [,4] [,5]
[1,] 9.660152e-01 5.557776e-07 2.474247e-01 -0.0185415248 8.139531e-01
[2,] 2.989917e-10 7.047560e-01 -1.357660e-09 0.0001775546 5.174325e-08
[3,] -3.734254e-01 3.671804e-06 2.608073e-01 -0.0990388112 3.230215e+00
[4,] 5.367914e-06 -6.728619e-05 -1.928126e-05 0.6273073468 5.587228e-04
[5,] 2.989064e-01 3.670783e-05 -7.924384e-01 -0.7506920943 1.689479e+01
Comparison between the linearized ODE and nonlinear ODE
(precise_coeff_5sp[[5]] - direct_coeff_5sp[[5]])/direct_coeff_5sp[[5]]
[,1] [,2] [,3] [,4] [,5]
[1,] -0.0004007798 7.320764e-03 0.0008132842 0.0011320063 0.0004001401
[2,] -0.0001265475 -3.913425e-08 0.0009166081 -0.0000665089 0.0001469353
[3,] -0.0004155999 4.155022e-03 0.0006365237 0.0001959350 -0.0001795353
[4,] -0.0003679745 7.152273e-05 0.0007680847 0.0007143146 0.0007743651
[5,] 0.0010231664 4.140094e-03 0.0019844202 0.0021181147 0.0019875209
Applying MDR S-map to the 5sp coupled foodchain model
Read dataset
da.range <- 1:200 # Subsample for data analysis
out.sample <- T # T/F for out-of-sample forecast
nout <- 2 # number of out-of-sample
da.name <- '5p_model'
do <- readRDS("rk4_result_5sp_model.obj")
dot <- do[da.range, 1] # data time
do <- do[da.range, -1] # subset of data without t
ndo <- nrow(do)
nin <- ndo - nout # library sample size
Standardization
# In-sample
do.mean <- apply(do[1:nin, ], 2, mean, na.rm = T) # mean abundance in in-sample
do.sd <- apply(do[1:nin, ], 2, sd, na.rm = T) # SD of abundance in in-sample
d <- do[1:(nin-1), ] # In-sample dataset at time t
d_tp1 <- do[2:(nin), ] # In-sample dataset at time t+1
ds <- (d - repmat(do.mean, nrow(d), 1))*repmat(do.sd, nrow(d), 1) ^ -1 # Normalized in-sample dataset at time t
ds_tp1 <- (d_tp1 - repmat(do.mean, nrow(d_tp1), 1))*repmat(do.sd, nrow(d_tp1), 1) ^ -1 # Normalized in-sample dataset at time t+1
# Out-sample
if(out.sample & nout != 0){
d.test <- do[nin:(ndo - 1), ] # Out-of-sample dataset at time t
dt_tp1 <- do[(nin + 1): ndo, ] # Out-of-sample dataset at time t+1
ds.test <- (d.test - repmat(do.mean, nrow(d.test), 1))*repmat(do.sd, nrow(d.test), 1) ^ -1 # Normalized out-of-sample dataset at time t
dst_tp1 <- (dt_tp1 - repmat(do.mean, nrow(dt_tp1), 1))*repmat(do.sd, nrow(dt_tp1), 1) ^ -1 # Normalized out-of-sample dataset at time t+1
}else{d.test <- dt_tp1 <- ds.test <- NULL}
# Compiled data at time t
ds.all <- rbind(ds, ds.test)
head(ds.all)
Initialization of pseudo-random numnber generator
seed <- 49563
set.seed(seed)
Find the optimal embedding dimension, based on univariate simplex projection
Emax <- 5 # equal to or smaller than the number of nodes (dimension of the ODE system)
cri <- 'rmse' # model selection
Ed <- NULL
forecast_skill_simplex <- NULL
for(i in 1:ncol(ds)){
spx.i <- simplex(ds[, i], E = 2:Emax)
Ed <- c(Ed, spx.i[which.min(spx.i[, cri])[1], 'E'])
forecast_skill_simplex <- c(forecast_skill_simplex, spx.i[which.min(spx.i[, cri])[1], 'rho'])
}
Ed # The optimal embedding dimension for each variable
[1] 2 3 3 3 2
Optimal embedding dimension and forecasting skill
Ed
[1] 2 3 3 3 2
forecast_skill_simplex
[1] 0.9597235 0.9789125 0.9540209 0.9699363 0.9646669
Find causal variables by CCM analysis that will be used for multiview embedding.
Warning: It is time consuming for calculating the causation for each node
CCM causality test for all node pairs (the function was updated at 2024/05/07)
ccm.out <- ccm.fast.demo(ds, Epair = T, cri = cri, Emax = Emax)
variable 1 ccm completed: 4.655 sec
variable 2 ccm completed: 5.451 sec
variable 3 ccm completed: 7.131 sec
variable 4 ccm completed: 5.187 sec
variable 5 ccm completed: 8.395 sec
ccm.sig <- ccm.out[['ccm.sig']]
ccm.rho <- ccm.out[['ccm.rho']]
#To avoid overwrite the original files, we save them with different names,identified by date.
DAY = Sys.Date()
saveRDS(ccm.sig, paste('ccm_sig', da.name,'nin', nin, DAY, '.obj', sep = '_'))
saveRDS(ccm.rho, paste('ccm_rho', da.name,'nin', nin, DAY, 'demo.obj', sep = '_'))
Just in case not conducting CCM but just loading the result of CCM at 2024/04/27
DAY = "2024-04-27"
#ccm.sig_old <- readRDS(paste('ccm_sig', da.name,'nin', nin, DAY, '.obj', sep = '_'))
#ccm.rho_old <- readRDS(paste('ccm_rho', da.name,'nin', nin, DAY, 'demo.obj', sep = '_'))
ccm.sig
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 1 1
[2,] 1 1 1 1 1
[3,] 1 0 1 1 1
[4,] 1 1 1 1 1
[5,] 1 1 1 1 1
#ccm.sig_old
ccm.rho
[,1] [,2] [,3] [,4] [,5]
[1,] 0.9926204 0.52102272 0.8795035 0.5119047 0.7022266
[2,] 0.4239670 0.99615225 0.6665958 0.8775080 0.6781655
[3,] 0.3827423 0.05067993 0.9962430 0.4746877 0.7138164
[4,] 0.6730854 0.49475755 0.5352258 0.9986359 0.7918114
[5,] 0.7003931 0.24789633 0.5921757 0.6243611 0.9988016
#ccm.rho_old
Based on the causal variables detected in CCM,
generates the time lag series (with maximum lag = 3),
conduct multivariate simplex projection and evaluate the forecast skill, and
selected the best kn(=100) SSR.
Perform multiview embedding analysis for each node
Warning: It is time consuming for running multiview embedding for each nodes
esele_lag <- esim.lag.demo(ds, ccm.rho, ccm.sig, Ed, kmax = 10000, kn = 100, max_lag = 3, Emax = Emax)
# To avoid overwrite the original files, we save them with different names, 'XXX_NEW'.
DAY = Sys.Date()
saveRDS(esele_lag, paste('eseleLag', da.name, 'nin', nin, DAY, 'demo.obj', sep='_'))
The list of the selected best multivariate SSRs for each node (species, sp) with
the correlation coefficient (rho)
the list of variables (X_1, X_2, …) used in each SSR, noting that the number of variables is identical to the optimal embedding dimension obtained in [6.2], and
the list of time lags (vlag_1, vlag_2, …) used in the corresponding variables (X_1, X_2, …).
as.data.frame(esele_lag)
dmatrix.mv <- mvdist.demo(ds, ds.all, esele_lag)
dmatrix.train.mvx <- dmatrix.mv[['dmatrix.train.mvx']]
dmatrix.test.mvx <- dmatrix.mv[['dmatrix.test.mvx']]
Leave-one-out cross-validation for finding the optimal parameters for MDR S-map analysis
–Warning: The cross-validation is the most time-consuming step in MDR S-map requiring massive computations and . Thus, we recommend dividing job into smaller parts (sub.da>1) or used parallel computation (parall=T, ncore>=1)
do.MDR.CV <- T
### The parameter cv.unit determines the precision of selected parameters and strongly influences computation time.
### In our cases, we used cv.unit=0.025 to obtain more precise estimations
### This parameter may be adjusted to 0.05 or even 0.1, depending on how sensitive the results to parameter precision.
cv.unit <- 0.025
alpha.so <- seq(0, 1, cv.unit); # Sequence of alpha
sub.da <- 1 # Divide the computation job into five parts
afsp <- eqsplit(1:length(alpha.so), sub.da) # Divide the parameter space based on alpha parameter
alf <- 1 # Run CV in the first parameter subset
# Cross-validation of MDR analysis
if(do.MDR.CV){
alpha.s <- alpha.so[afsp[alf, 1]:afsp[alf, 2]] # Subset parameter space
cv.ind <- cv.MDR.demo(ds, ds_tp1, dmatrix.list = dmatrix.train.mvx,
parall = T, ncore = 24, keep_intra = T, alpha.seq = alpha.s)
# To avoid overwrite the original files, we save them with different names, 'XXX_NEW'.
saveRDS(cv.ind, paste(da.name, 'nin', nin, 'cvunit', cv.unit, 'alph', alpha.s[1]*100, DAY, 'cvout_Nmvx_Rallx.obj', sep = '_'))
}
Select the optimal parameter set for regularized regression model with the minimal MSE
paracv.demo <- secv.demo(cv.ind)
paracv.demo
Fitting MDR S-map based on the parameters selected by CV
Note that the linear regression model in MDR S-map includes all of variables as X even if some of them were not selected as the causal variables at the step of using CCM. This is reasonable because S-map coefficients represent the effects of variables at short timescale while CCM evaluates the causality in the whole period of the time series.
#Setting
do.MDR <- F
cv.unit <- 0.025
ptype <- 'aenet' #enet:elastic-net or msaenet: adaptive elastic-net
#Fitting the MDR S-map
smap.5sp <- MDRsmap.demo(paracv = paracv.demo, ptype = ptype, keep_intra = T, out.sample = T,
ds,ds_tp1,ds.test,dst_tp1,
dmatrix.list = dmatrix.train.mvx,
dmatrix.test.list = dmatrix.test.mvx)
In case for saving the results
DAY = Sys.Date()
nr.out <- smap.5sp[['nr.out']];
saveRDS(nr.out, paste(da.name,'_nin', nin, 'cvunit', cv.unit, ptype, DAY, 'nrout_Nmvx_Rallx_demo_NEW.obj',sep = '_'))
# Save interaction Jacobian matrices at all time points
saveRDS(smap.5sp[['jcof']], paste(da.name,'nin', nin, 'cvunit', cv.unit, ptype, DAY, '_jcof_Nmvx_Rallx_demo_NEW.obj', sep='_'))
S-map coefficients
smap.5sp[['jcof']]
#or
#"5p_model_nin_198_cvunit_0.025_aenet_2024-07-19__jcof_Nmvx_Rallx_demo_NEW.obj"
Loading data
smap.5sp <- list()
smap.5sp[['jcof']] <- readRDS("5p_model_nin_198_cvunit_0.025_aenet_2024-07-19__jcof_Nmvx_Rallx_demo_NEW.obj")
smap.5sp[['jcof']]
CRIS (Eqn.21) and direct_IS (Eqn.19)
#normalization by SD
sd_adj_35 <- as.numeric(rk4_result.sd[5] / rk4_result.sd[3])
CRIS_dC1dC1 <- precise_coeff_5sp[[1]][3,3] - 1.0
direct_dC1dC1 <- direct_coeff_5sp[[1]][3,3] - 1.0
CRIS_dC1dR <- sd_adj_35 * precise_coeff_5sp[[1]][3,5]
direct_dC1dR <- sd_adj_35 * direct_coeff_5sp[[1]][3,5]
for(j in 2:200) {
CRIS_dC1dC1 <- append(CRIS_dC1dC1, precise_coeff_5sp[[j]][3,3] - 1.0)
direct_dC1dC1 <- append(direct_dC1dC1, direct_coeff_5sp[[j]][3,3] - 1.0)
CRIS_dC1dR <- append(CRIS_dC1dR, sd_adj_35 * precise_coeff_5sp[[j]][3,5])
direct_dC1dR <- append(direct_dC1dR, sd_adj_35 * direct_coeff_5sp[[j]][3,5])
}
Dataframe for basic plot
final_ID <- 201
tau <- 5.0
food_chain_result <- data.frame(time = rk4_result$t[-final_ID], std_smapC1toC1 = rk4_m_smap_res2_s$smap_coefficients[[1]]$c_3[-final_ID] - 1.0, std_smapRtoC1 = rk4_m_smap_res2_s$smap_coefficients[[1]]$c_5[-final_ID], mdr_smapC1toC1 = subset(smap.5sp[['jcof']], variable == j_C1)[1:200, 7] - 1.0, mdr_smapRtoC1 = subset(smap.5sp[['jcof']], variable == j_C1)[1:200, 9], IIS_dC1dC1 = tau*dC1dC1(rk4_result)[-final_ID], IIS_dC1dR = tau*dC1dR(rk4_result)[-final_ID], CRIS_dC1dC1, CRIS_dC1dR, direct_dC1dC1, direct_dC1dR)
food_chain_result
Dataframe for ggplot
#For lineplots
rk4_C1toC1_ggplot <- data.frame(
time = food_chain_result$time,
IntS = food_chain_result$std_smapC1toC1,
type = rep("01std_smapC1toC1", 200)
)
rk4_C1toC1_ggplot <- rbind.data.frame(
rk4_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$mdr_smapC1toC1,
type = rep("02mdr_smapC1toC1", 200)
)
)
rk4_C1toC1_ggplot <- rbind.data.frame(
rk4_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$CRIS_dC1dC1,
type = rep("03CRIS_dC1dC1", 200)
)
)
rk4_C1toC1_ggplot <- rbind.data.frame(
rk4_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$IIS_dC1dC1,
type = rep("04IIS_dC1dC1", 200)
)
)
rk4_RtoC1_ggplot <- data.frame(
time = food_chain_result$time,
IntS = food_chain_result$std_smapRtoC1,
type = rep("01std_smapRtoC1", 200)
)
rk4_RtoC1_ggplot <- rbind.data.frame(
rk4_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$mdr_smapRtoC1,
type = rep("02mdr_smapRtoC1", 200)
)
)
rk4_RtoC1_ggplot <- rbind.data.frame(
rk4_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$CRIS_dC1dR,
type = rep("03CRIS_dC1dR", 200)
)
)
rk4_RtoC1_ggplot <- rbind.data.frame(
rk4_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$IIS_dC1dR,
type = rep("04IIS_dC1dR", 200)
)
)
#For error plots
rk4_error_C1toC1_ggplot <- data.frame(
time = food_chain_result$time,
IntS = food_chain_result$std_smapC1toC1 - food_chain_result$direct_dC1dC1,
type = rep("01std_smapC1toC1", 200)
)
rk4_error_C1toC1_ggplot <- rbind.data.frame(
rk4_error_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$mdr_smapC1toC1 - food_chain_result$direct_dC1dC1,
type = rep("02mdr_smapC1toC1", 200)
)
)
rk4_error_C1toC1_ggplot <- rbind.data.frame(
rk4_error_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$CRIS_dC1dC1 - food_chain_result$direct_dC1dC1,
type = rep("03CRIS_dC1dC1", 200)
)
)
rk4_error_C1toC1_ggplot <- rbind.data.frame(
rk4_error_C1toC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$IIS_dC1dC1 - food_chain_result$direct_dC1dC1,
type = rep("04IIS_dC1dC1", 200)
)
)
rk4_error_RtoC1_ggplot <- data.frame(
time = food_chain_result$time,
IntS = food_chain_result$std_smapRtoC1 - food_chain_result$direct_dC1dR,
type = rep("01std_smapRtoC1", 200)
)
rk4_error_RtoC1_ggplot <- rbind.data.frame(
rk4_error_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$mdr_smapRtoC1 - food_chain_result$direct_dC1dR,
type = rep("02mdr_smapRtoC1", 200)
)
)
rk4_error_RtoC1_ggplot <- rbind.data.frame(
rk4_error_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$CRIS_dC1dR - food_chain_result$direct_dC1dR,
type = rep("03CRIS_dC1dR", 200)
)
)
rk4_error_RtoC1_ggplot <- rbind.data.frame(
rk4_error_RtoC1_ggplot, data.frame(
time = food_chain_result$time,
IntS = food_chain_result$IIS_dC1dR - food_chain_result$direct_dC1dR,
type = rep("04IIS_dC1dR", 200)
)
)
plot(food_chain_result$time, food_chain_result$CRIS_dC1dC1, type = "l", col = "#f6aa00", xlim = c(1400, 1800), ylim = c(-3, 2.0), xlab = "time", ylab = "Smap coefficient: C1 -> C1")
par(new = T)
plot(food_chain_result$time, food_chain_result$std_smapC1toC1, type = "l", col = "#ff4b00", xlim = c(1400, 1800), ylim = c(-3, 2.0), xlab = "", ylab = "")
par(new = T)
plot(food_chain_result$time, food_chain_result$mdr_smapC1toC1, type = "l", col = "#4dc4ff", xlim = c(1400, 1800), ylim = c(-3, 2.0), xlab = "", ylab = "")
par(new = T)
plot(food_chain_result$time, food_chain_result$IIS_dC1dC1, type = "l", lty = "dashed", col = "#804000", xlim = c(1400, 1800), ylim = c(-3, 2.0), xlab = "", ylab = "")
fig3a_ggplot <- ggplot(data = rk4_C1toC1_ggplot) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00", "#804000")) +
scale_shape_manual(values = c(15, 16, 17, 5)) +
 geom_line(aes(x = time, y = IntS, group = type, color = type)) +
geom_point(aes(x = time, y = IntS, color = type, shape = type), size = 2) +
xlim(1400, 1800) + theme_bw() + theme(legend.position = "bottom") +
labs(title = "Time evolution: C1 -> C1", x = "time", y = "C1 -> C1")
fig3b_ggplot <- ggplot(rk4_error_C1toC1_ggplot, aes(type, abs(IntS))) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00", "#804000")) +
scale_shape_manual(values = c(15, 16, 17, 5)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(shape = type, color = type, alpha = 0.1), width = 0.25, size = 3) +
theme_bw() + theme(legend.position = "bottom") +
labs (title = "Comparison to direct_evaluation: C1 -> C1")
plot(fig3a_ggplot)
ggsave("fig3a.pdf", width = 6, height = 4)
plot(fig3b_ggplot)
ggsave("fig3b.pdf", width = 4, height = 4)
Difference between directly-evaluated CIS and CIS by linearized ODE
summary(abs(food_chain_result$CRIS_dC1dC1 - food_chain_result$direct_dC1dC1))
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.157e-06 7.075e-05 2.301e-04 3.547e-04 5.008e-04 1.844e-03
fig3c_ggplot <- ggplot(data = rk4_RtoC1_ggplot) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00", "#804000")) +
scale_shape_manual(values = c(15, 16, 17, 5)) +
 geom_line(aes(x = time, y = IntS, group = type, color = type)) +
geom_point(aes(x = time, y = IntS, color = type, shape = type), size = 2) +
xlim(1400, 1800) + theme_bw() + theme(legend.position = "bottom") +
labs(title = "Time evolution: R -> C1", x = "time", y = "R -> C1")
fig3d_ggplot <- ggplot(rk4_error_RtoC1_ggplot, aes(type, abs(IntS))) +
scale_colour_manual(values = c("#ff4b00", "#4dc4ff", "#f6aa00", "#804000")) +
scale_shape_manual(values = c(15, 16, 17, 5)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(aes(shape = type, color = type, alpha = 0.1), width = 0.25, size = 3) +
theme_bw() + theme(legend.position = "bottom") +
labs (title = "Comparison to direct_evaluation: R -> C1")
plot(fig3c_ggplot)
ggsave("fig3c.pdf", width = 6, height = 4)
plot(fig3d_ggplot)
ggsave("fig3d.pdf", width = 4, height = 4)
Difference between directly-evaluated CIS and CIS by linearized ODE
summary(abs(food_chain_result$CRIS_dC1dR - food_chain_result$direct_dC1dR))
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.000e-09 5.918e-06 5.095e-05 5.600e-04 3.540e-04 1.091e-02
#for diagonal element (C1 -> C1)
plot(food_chain_result$std_smapC1toC1 ~ food_chain_result$direct_dC1dC1,
xlab = "direct evaluation (C1 -> C1)",
ylab = "standard S-map coefficient (C1 -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$std_smapC1toC1[-200], food_chain_result$direct_dC1dC1[-200])
[1] 0.8374068
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$std_smapC1toC1[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$std_smapC1toC1[-200])
Residuals:
Min 1Q Median 3Q Max
-1.05026 -0.17010 -0.00113 0.16673 1.45238
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02650 0.02399 -1.105 0.271
food_chain_result$std_smapC1toC1[-200] 0.91642 0.04262 21.504 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3337 on 197 degrees of freedom
Multiple R-squared: 0.7013, Adjusted R-squared: 0.6997
F-statistic: 462.4 on 1 and 197 DF, p-value: < 2.2e-16
plot(food_chain_result$mdr_smapC1toC1 ~ food_chain_result$direct_dC1dC1,
xlab = "direct evaluation (C1 -> C1)",
ylab = "MDR S-map coefficient (C1 -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$mdr_smapC1toC1[-200], food_chain_result$direct_dC1dC1[-200])
[1] 0.7521856
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$mdr_smapC1toC1[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$mdr_smapC1toC1[-200])
Residuals:
Min 1Q Median 3Q Max
-1.27802 -0.08214 0.05370 0.19339 0.94716
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.12260 0.03070 -3.994 9.17e-05 ***
food_chain_result$mdr_smapC1toC1[-200] 1.01376 0.06327 16.022 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4023 on 197 degrees of freedom
Multiple R-squared: 0.5658, Adjusted R-squared: 0.5636
F-statistic: 256.7 on 1 and 197 DF, p-value: < 2.2e-16
plot(food_chain_result$IIS_dC1dC1 ~ food_chain_result$direct_dC1dC1,
xlab = "direct evaluation (C1 -> C1)",
ylab = "instantaneous interaction strength (C1 -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$IIS_dC1dC1[-200], food_chain_result$direct_dC1dC1[-200])
[1] 0.7062711
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$IIS_dC1dC1[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dC1[-200] ~ food_chain_result$IIS_dC1dC1[-200])
Residuals:
Min 1Q Median 3Q Max
-1.67035 -0.13355 0.06334 0.30644 1.02468
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.005217 0.030881 0.169 0.866
food_chain_result$IIS_dC1dC1[-200] 0.729854 0.052123 14.003 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4322 on 197 degrees of freedom
Multiple R-squared: 0.4988, Adjusted R-squared: 0.4963
F-statistic: 196.1 on 1 and 197 DF, p-value: < 2.2e-16
#for off-diagonal element (R -> C1)
plot(food_chain_result$std_smapRtoC1 ~ food_chain_result$direct_dC1dR,
xlab = "direct evaluation (R -> C1)",
ylab = "standard S-map coefficient (R -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$std_smapRtoC1[-200], food_chain_result$direct_dC1dR[-200])
[1] 0.5088147
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dR[-200] ~ food_chain_result$std_smapRtoC1[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dR[-200] ~ food_chain_result$std_smapRtoC1[-200])
Residuals:
Min 1Q Median 3Q Max
-1.3354 -0.2343 -0.1558 -0.0639 4.7867
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.11222 0.07033 1.596 0.112
food_chain_result$std_smapRtoC1[-200] 0.96883 0.11679 8.296 1.68e-14 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.7805 on 197 degrees of freedom
Multiple R-squared: 0.2589, Adjusted R-squared: 0.2551
F-statistic: 68.82 on 1 and 197 DF, p-value: 1.678e-14
plot(food_chain_result$mdr_smapRtoC1 ~ food_chain_result$direct_dC1dR,
xlab = "direct evaluation (R -> C1)",
ylab = "MDR S-map coefficient (R -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$mdr_smapRtoC1[-200], food_chain_result$direct_dC1dR[-200])
[1] 0.2223845
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dR[-200] ~ food_chain_result$mdr_smapRtoC1[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dR[-200] ~ food_chain_result$mdr_smapRtoC1[-200])
Residuals:
Min 1Q Median 3Q Max
-0.6196 -0.4075 -0.2765 -0.0205 5.1085
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.22700 0.09899 2.293 0.02290 *
food_chain_result$mdr_smapRtoC1[-200] 0.59385 0.18549 3.201 0.00159 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.884 on 197 degrees of freedom
Multiple R-squared: 0.04945, Adjusted R-squared: 0.04463
F-statistic: 10.25 on 1 and 197 DF, p-value: 0.001594
plot(food_chain_result$IIS_dC1dR ~ food_chain_result$direct_dC1dR,
xlab = "direct evaluation (R -> C1)",
ylab = "instantaneous interaction strength (R -> C1)",
)
abline(coef = c(0,1))
cor(food_chain_result$IIS_dC1dR[-200], food_chain_result$direct_dC1dR[-200])
[1] 0.7100644
#Y should be the theoretical coefficient
summary(lm(food_chain_result$direct_dC1dR[-200] ~ food_chain_result$IIS_dC1dR[-200]))
Call:
lm(formula = food_chain_result$direct_dC1dR[-200] ~ food_chain_result$IIS_dC1dR[-200])
Residuals:
Min 1Q Median 3Q Max
-2.3629 -0.1519 -0.1442 -0.0614 4.1714
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.15031 0.05065 2.967 0.00338 **
food_chain_result$IIS_dC1dR[-200] 0.43744 0.03091 14.154 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6384 on 197 degrees of freedom
Multiple R-squared: 0.5042, Adjusted R-squared: 0.5017
F-statistic: 200.3 on 1 and 197 DF, p-value: < 2.2e-16
Note that:
delta_33 <- 1
precise_data <- data.frame(t = rk4_result$t[1], precise_dC1dC1 = precise_coeff_5sp[[1]][3,3] - delta_33, inst_precise_dC1dC1 = inst_precise_coeff_5sp[[1]][3,3] - delta_33, tau0.5_precise_dC1dC1 = tau0.5_precise_coeff_5sp[[1]][3,3] - delta_33, tau1.0_precise_dC1dC1 = tau1.0_precise_coeff_5sp[[1]][3,3] - delta_33, tau2.0_precise_dC1dC1 = tau2.0_precise_coeff_5sp[[1]][3,3] - delta_33, tau10.0_precise_dC1dC1 = tau10.0_precise_coeff_5sp[[1]][3,3] - delta_33)
for(j in 2:200) precise_data <- rbind.data.frame(precise_data, c(rk4_result$t[j], precise_coeff_5sp[[j]][3,3] - delta_33, inst_precise_coeff_5sp[[j]][3,3] - delta_33, tau0.5_precise_coeff_5sp[[j]][3,3] - delta_33, tau1.0_precise_coeff_5sp[[j]][3,3] - delta_33, tau2.0_precise_coeff_5sp[[j]][3,3] - delta_33, tau10.0_precise_coeff_5sp[[j]][3,3] - delta_33))
plot(rk4_result$t[101:181], tau2*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.005, 0.005), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 0.01", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$inst_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.005, 0.005), xlab = "", ylab = "")
legend("bottomright",
legend = c("instantaneous interaction strength", "cumulative interaction strength"),
col = c("black", "#ff4b00"),
lty = c(1, 2),
cex = c(0.8)
)
plot(rk4_result$t[101:181], tau3*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.4, 0.3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 0.5", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau0.5_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.4, 0.3), xlab = "", ylab = "")
plot(rk4_result$t[101:181], tau4*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(0.0, 0.15), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 1.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau1.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(0.0, 0.15), xlab = "", ylab = "")
plot(rk4_result$t[101:181], tau5*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.4, 0.3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 2.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau2.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.4, 0.3), xlab = "", ylab = "")
plot(rk4_result$t[101:181], tau1*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-2, 2), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 5.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-2, 2), xlab = "", ylab = "")
plot(rk4_result$t[101:181], tau6*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-2, 3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 10.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau10.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-2, 3), xlab = "", ylab = "")
pdf("fig4a.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau2*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.005, 0.005), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 0.01", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$inst_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.005, 0.005), xlab = "", ylab = "")
legend("bottomright",
legend = c("instantaneous interaction strength", "cumulative interaction strength"),
col = c("black", "#ff4b00"),
lty = c(1, 2),
cex = c(0.8)
)
dev.off()
null device
1
pdf("fig4b.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau3*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.4, 0.3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 0.5", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau0.5_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.4, 0.3), xlab = "", ylab = "")
dev.off()
null device
1
pdf("fig4c.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau4*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(0.0, 0.15), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 1.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau1.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(0.0, 0.15), xlab = "", ylab = "")
dev.off()
null device
1
pdf("fig4d.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau5*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-0.4, 0.3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 2.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau2.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-0.4, 0.3), xlab = "", ylab = "")
dev.off()
null device
1
pdf("fig4e.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau1*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-2, 2), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 5.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-2, 2), xlab = "", ylab = "")
dev.off()
null device
1
pdf("fig4f.pdf", width = 4, height = 4)
plot(rk4_result$t[101:181], tau6*dC1dC1(rk4_result)[101:181], type = "l", col = "black", ylim = c(-2, 3), xlab = "time", ylab = "Interaction strengths C1 -> C1", lwd = 1, main = expression(paste(tau, "= 10.0", sep="")))
par(new = T)
plot(rk4_result$t[101:181], precise_data$tau10.0_precise_dC1dC1[101:181], type = "l", col = "#ff4b00", lwd = 2, lty = 2, ylim = c(-2, 3), xlab = "", ylab = "")
dev.off()
null device
1